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SECTION  I 
INTRODUCTION 


This  first  annual  technical  report  documents  the  current 
status  and  accomplishments  of  the  research  performed  by  the 
University  of  Dayton  for  the  Air  Force  Office  of  Scientific 
Research  (AFOSR) ,  Fast  Algorithms  Initiative.  The  research 
documented  herein  was  carried  out  by  the  University  of  Dayton 
Research  Institute  (UDRI)  as  the  prime  contractor  and  by  the 
University  of  Pittsburgh  Institute  for  Computational  Mathematics 
and  Applications  (ICMA)  as  the  subcontractor. 

1 . 1  BACKGROUND 

The  so-called  turbulence  problem  in  fluid  mechanics  remains 
unsolved,  but  not  for  a  lack  of  theoretical  Interest  or  of 
practical  Importance.  From  the  theoretical  viewpoint,  the 
understanding  and  prediction  of  turbulent  shear  flows  from  first 
principles  have  remained  elusive  and  continue  to  provide 
challenging  problems  in  fluid  dynamics,  computational  physics  and 
mathematics,  and  numerical  algorithm  development.  From  a 
practical  viewpoint,  many  research  and  development  problems  of 
Interest  to  the  Air  Force,  ranging  from  combustion  in  jet  engines 
and  rocket  motors  to  aerodynamics  of  wings  and  surfaces,  give 
rise  to  the  transport  of  mass,  momentum,  and  energy  that  are 
predominantly  governed  by  turbulent  flows,  thereby  necessitating 
their  understanding  and  prediction. 

It  is  generally  accepted  that  fluid  mechanical  turbulence 
begins  in  an  erstwhile  laminar  flowfleld  as  the  onset  and  growth 
of  hydrodynamic  instabilities.  According  to  this  traditional 
viewpoint,  as  the  energy  driving  the  macroscopic  instabilities 
cascades  through  higher  and  higher  wavenumbers,  culminating  in 
what  is  commonly  referred  to  as  the  fully  developed  turbulence, 
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the  overall  flowfleld  Is  characterized  by  the  interaction  of  a 
very  large  number  of  disparate  length  and  time  scales.  This 
viewpoint  has  led  to  the  predictive  approach  which  seeks  the 
evolution  of  turbulence  from  the  solution  of  the  time-dependent 
Navler-Stokes  equations.  If  these  equations  can  be  solved 
exactly  for  the  entire  range  of  the  time  and  length  scales 
prevailing  In  a  realistic  flowfleld,  then  the  resulting  solutions 
explicitly  describe  and  predict  the  turbulent  flow. 

Unfortunately,  this  approach  has  been  thwarted  by  the  enormous 
computing  requirements  imposed  by  the  time-dependent, 
multidimensional  Navler-Stokes  equations  In  simulating  a 
realistic  flow,  which  requirements  are  well  beyond  the 
capabilities  of  today's  largest  and  fastest  computers. 

The  only  feasible  approach  for  turbulent  flow  predictions  at 
present  comprises  the  direct  computation  of  the  complete 
equations  on  a  coarsely  resolved  grid  (dictated  by  the  available 
computing  resources)  for  describing  the  "large "-scale  motions, 
the  development  of  model (s)  for  the  "small"  subgrid-scale  (SGS) 
turbulent  mixing,  and  the  proper  coupling  of  the  SGS  modeling  to 
the  large-eddy  simulation  (LBS).  Indeed,  the  prospects  for 
economical  and  accurate  LES  of  turbulent  shear  flows  have  been 
greatly  enhanced  with  the  Increasing  availability  of  vector 
computers  of  greater  and  greater  size  and  speed,  the  potential 
for  major  Improvements  in  the  development  of  efficient  numerical 
algorithms,  and  the  considerable  scope  for  efficient  matching  of 
these  algorithms  to  emerging  computer  architectures. 

1.2  SCOPE  AND  OBJECTIVES 

The  overall  goal  of  this  UDRI-ICMA  joint  research  Is  the 
computational  Investigation  of  the  turbulent  mixing  layers  of 
several  prototype  circular- jet  configurations.  Initial  attention 
is  directed  at  an  examination  of  a  single  free  jet  expanding  Into 
a  quiescent  environment,  with  emphasis  on  the  development  of 


[2] 


suitable  SGS  turbulence  models  and  the  computation  of  large-scale 
motions.  Specific  areas  addressed  during  the  first  year  Include 
the  following: 

•  Investigation  of  the  axlsymmetrlc  mixing  layer  of 
circular  jets  with  LES  and  SGS  modeling. 

•  Investigation  and  analysis  of  variable  reduction 
algorithms,  such  as  the  dual  variable  method. 

•  Investigation  and  analysis  of  numerical  algorithms  based 
on  weighted  combinations  of  stable  and  unstable  schemes. 

•  Investigation  of  vectorlzatlon  of  solution  algorithms  and 
its  efficient  Implementation. 

Our  research  to  date  has  entailed  a  critical  review  and 
assessment  of  the  LES  and  SGS  modeling  procedures  prior  to  the 
selection  of  the  specific  avenues  that  have  been  explored. 
Furthermore,  the  research  has  Involved  several  Interdependent 
phases  consisting  of  a  number  of  tasks.  Some  of  these  tasks  are 
independently  conducted  by  UDRI  (e.g.,  SGS  turbulence  modeling) 
and  by  ICMA  (e.g.,  modification  and  adaptation  of  an  existing 
two-dimensional  procedure)  and  some  have  entailed  joint  efforts 
(e.g.,  boundary-condition  modifications  and  algebraic-turbulence- 
model  incorporation  and  testing) .  These  are  discussed  In  the 
next  section. 

1.3  OUTLINE  OF  REPORT 


Section  II  discusses  our  accomplishments  to  date,  and  gives 
the  current  status  of  our  research.  Section  III  presents  the 
documentation  from  research  sponsored  under  this  program.  The 
research  personnel  at  UDRI  and  ICMA  supported  by  this  program  are 
shown  in  Section  IV. 
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SECTION  II 
STATUS  07  RESEARCH 


There  has  been  significant  progress  in  achieving  the 
research  objectives  outlined  in  the  previous  section.  A  brief 
summary  of  our  accomplishments  to  date  is  presented  below. 

2.1  MATHEMATICAL  DESCRIPTION 

The  prediction  of  jet  turbulence  in  the  present  research  has 
adopted  the  traditional  viewpoint  that  it  is  based  upon  the 
solutions  of  the  time-dependent  Navier-Stokes  equations.  It  must 
be  emphasized,  nevertheless,  that  these  equations  remain  to  be 
fully  tested  for  completeness  by  computations.  Be  that  as  it 
may,  the  time-averaged  information  in  a  turbulent  jet  is  obtained 
by  solving  the  time-dependent  equations  and  averaging  the 
solutions.  The  prediction  of  the  large-scale  motion  by  LES 
requires  the  computation  of  the  three-dimensional  equations, 
since  the  large  eddies  are  highly  anisotropic.  Thus,  while  the 
circular-jet  geometry  and  its  mean-flow  characteristics  preserve 
axlsymmetry,  such  an  assumption  for  LES  is  invalid.  Within  the 
constraints  of  the  overall  research  program,  however,  the  present 
predictive  research  is  based  upon  the  two-dimensional 
formulation.  As  regards  the  small-scale  motion,  even  with  our 
calculation  methodology  of  LES  the  smallest  scale  explicitly 
resolved  is  very  much  larger  than  the  Kolmogoroff  [1940]  scale 
for  viscous  dissipation.  Thus,  while  the  smallest  scales  are 
truly  Isotropic  and  possess  universal  characteristics,  that  is 
not  the  case  at  the  (arbitrary)  LES  cutoff  length  scale.  This 
necessitates  a  closure  model  for  wavenumbers  larger  than  that 
corresponding  to  the  LES  cutoff,  with  the  consequent 
indeterminacy  and  the  ad  hoc  introduction  of  turbulence  modeling. 
The  present  SGS  modeling  approach  is  crude  but  provides  a  useful 
starting  point  for  the  incorporation  of  a  variable  viscosity  in 
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our  calculation  procedure.  These  LES  and  SGS  turbulence  modeling 
approaches  are  discussed  next. 

2.1.1  Large-Scale  Motion 

The  Initial  research  effort  has  addressed  the  assessment  and 
analysis  of  LES  procedures  with  a  view  to  arriving  at  suitable 
schemes  for  investigating  circular- jet  turbulence. 

The  basis  for  our  predictive  computation  is  the  viewpoint 
that  turbulence  consists  of  chaotic  motion  (with  persistent 
organized  structures)  at  a  range  of  length  scales  that  increases 
rapidly  with  the  Reynolds  number.  This  range  of  scales  stems 
from  the  convective  terms  of  the  Navler-Stokes  equations.  The 
entire  range  of  scales  can  be  numerically  resolved  with  no 
modification  of  the  convective  terms  only  for  sufficiently  low 
Reynolds  numbers  (at  which  complete  numerical  simulation  is 
possible  with  currently  available  computers,  but  turbulence  can 
not  be  maintained) .  When  the  desired  scale  range  exceeds  the 
lower  bound  imposed  by  computer  capacity  (as  in  LES),  the  scales 
that  are  smaller  than  the  LES  cutoff  are  discarded,  and  the 
influence  of  these  discarded  scales  upon  the  retained  scales  must 
be  modeled.  This  necessitates  the  modification  of  the  governing 
equations,  since  the  computations  involving  the  unmodified 
Navler-Stokes  equations  and  the  arbitrary  cutoff  in  the 
resolution  have  no  relation  to  the  real  fluid  physics.  This  lack 
of  realism  can  manifest  itself  in  several  ways: 

•  The  numerical  algorithm  for  the  unmodified  LES  equations 
can  become  unstable  as  the  smallest  resolved  scales 
accumulate  energy. 

•  The  use  of  energy-conserving  numerical  approximations  can 
lead  to  a  nonphysical  equilibrium  distribution  of  the 
energy  among  the  finite  degrees  of  freedom. 
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•  When  the  viscosity  is  not  zero  but  is  too  small  to  allow 
accurate  resolutions  of  the  dissipation  scales,  an 
energy-conserving  algorithm  collects  energy  at  the 
smallest  resolved  scales  until  there  is  equilibrium 
between  the  dissipation  and  cascade  rates.  Indeed,  it 
has  been  demonstrated  that  this  excess  energy,  trapped  at 
the  mesh  scale  rather  than  cascading  to  the  Kolmogoroff 
scale,  produces  too  rapid  an  energy  transfer  from  the 
large  scales.  Such  an  occurrence  can  well  be  expected  if 
the  small  scales  act  on  the  large  scales  as  an  eddy 
viscosity  with  a  value  (proportional  to  the  length  and 
velocity  scales  of  the  trapped  energy)  that  is  enhanced 
by  the  entrapment. 

Thus,  one  of  the  most  important  modifications  of  the  Navier- 
Stokes  equations  prior  to  current  LES  computation  is  the 
Incorporation  of  terms  that  provide  the  mechanism  for  removal  of 
energy  from  the  computed  scales  that  mimics  as  closely  as 
possible  the  physical  cascade  process.  This  aspect  is  currently 
under  investigation. 

Concurrently,  LES  investigations  of  an  approximate  nature 
have  been  initiated  for  addressing  the  single  free  jet  issuing 
into  an  unbounded  ambient  domain.  This  research  is  based  upon 
the  ALgorithms  for  GAs  Equations  (ALGAE),  which  is  a  general- 
purpose  computational  procedure  developed  at  ICMA.  The  version 
of  ALGAE  available  at  the  start  of  this  research  program  is 
capable  of  simulating  transient.  Incompressible  or  compressible, 
laminar  flows  in  two-dimensional  planar  and  axisymmetric 
geometries.  ALGAE  makes  use  of  the  dual  variable  method  for 
reducing  the  number  of  difference  equations  that  need  to  be 
solved  at  each  time  Interval.  The  availability  of  this 
particular  feature  and  the  potential  for  extension  to  address 
three-dimensional  flows  in  future  have  led  to  our  choice  of  ALGAE 
for  Investigating  the  large-scale  motions,  despite  our 
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recognition  of  its  current  limitations  in  providing  true  LES  of 
the  circular- Jet  turbulence.  Our  present  intent  has  been  to 
treat  ALGAE  as  a  baseline  procedure  that  can  be  refined  and 
optimized  for  the  problem  at  hand  through  successive 
implementation  of  a  number  of  modifications.  The  modifications 
addressed  during  this  research  are  discussed  subsequently. 

2.1.2  Small-Scale  Motion 

Several  essential  requirements  must  be  met  by  the  SGS 
turbulence  model,  if  it  is  to  provide  the  correct  asymptotic 
decoupling  from  the  LES  computations.  To  start  with,  SGS 
model (s)  must  satisfy  certain  matching  requirements  at  two 
distinguished  limits  of  wavenumbers.  These  arise  from  the  fact 
that  the  essential  physics  of  turbulence  is  continuously 
distributed  over  the  entire  spectrum  of  wavenumbers,  ranging  from 
k  *  0  to  k  *  (k^  is  the  dissipation  wavenumber  corresponding 

to  the  Kolmogoroff  microscale,  and  wavenumbers  exceeding  kd 
belong  to  the  realm  of  molecular  motion  and  do  not  concern  us) , 
whereas  the  division  of  wavenumbers  (or  length  scales)  into 
"large  eddies"  and  modeled  "subgrid  quantities”  is  simply 
arbitrary.  Thus,  the  SGS  model  prediction  must  merge  smoothly  to 
the  "laminar"  limit  at  k^  on  the  one  hand,  and  to  the  "grid-size" 
limit  (which  is  the  largest  resolvable  wavenumber  in  LES)  on 
the  other.  Unfortunately  k^  can,  and  often  does,  change  with 
available  computing  power  (and  also  when,  for  example,  the 
numerical  resolution  is  doubled  or  halved) .  Therefore  the  SGS 
model  predictions  must  also  remain  Invariant  of  these  macroscopic 
changes,  when  they  are  coupled  to  the  LES  computations. 

SGS  model (s)  must  go  beyond  the  mere  matching  requirements. 
It  is  essential  that  they  also  provide  for  the  phenomenological 
effects  of  the  laminar  viscosity,  dlffuslvlty,  etc.  (which 
represent  the  averages  over  the  molecular  dynamics) . 


Furthermore,  additional  coefficients  of  viscosity,  diffusivity, 
etc.  for  the  eddy  motion  are  seen  to  arise  from  the  averaging 
over  the  unresolved  wavenumbers  ( that  are  bounded  from  below  by 
and  from  above  by  k^) .  The  effects  of  thes^  coefficients  on 
the  large-scale  motion  must  be  incorporated  as  well  by  the  SGS 
model.  The  model (s)  must,  of  course,  predict  the  onset  of 
initial  instability  and  its  subsequent  growth  in  an  erstwhile 
laminar  flow,  with  an  increase  in  the  parameter  --  parameters  of 
interest  {e.g.,  Reynolds  number).  Finally,  the  model  must  be 
computationally  efficient  when  it  is  coupled  to  the  LES 
procedure . 

The  nature  of  small-scale  motion  could  be  further  elucidated 

by  examining  an  important  field  variable  of  turbulence.  This  is 

the  local  dissipation  rate  per  unit  mass  defined  as  e*  (v/2) 

2 

( 3u^/3Xj  +  3Uj/axi)  ,  where  e  is  averaged  over  volumes  of  0(L),  L 
being  the  energy  scale  of  motion,  with  L  >  1/k^;  x  and  u  denote 
the  coordinate  and  velocity,  and  v  *  p/p  is  the  kinematic 
viscosity,  with  p  and  p  respectively  denoting  the  mass  density 
and  molecular  viscosity  coefficient.  A  characteristic  property 
of  turbulence  is  that  e  increasingly  varies  with  increasing 
Reynolds  number.  This  random  character  of  e  was  suggested  by 
Kolmogoroff  [1962]  and  Obukhoff  [1962]  in  a  refinement  of 
Kolmogorof f ' s  original  hypotheses  [1940].  They  argued  that  e 
should  not  only  be  random  but  that  its  distribution  be  lognormal. 
Yaglom  [1966]  demonstrated  that  the  assumption  of  a  cascade 
process  of  energy  transfer  from  the  very  large-  to  the  very 
small-scale  motions  implies  that  ln(e)  should  be  Gaussian  if  the 
transfer  stages  are  statistically  similar  and  Independent.  This 
lognormality  prediction  by  Kolmogoroff,  Obukhoff,  and  Yaglom  has 
been  subsequently  verified  by  a  number  of  experimental 
investigations.  Thus,  while  the  large-scale  motions  are  properly 
treated  deterministically  through  LES,  the  SGS  turbulence 
model (s)  must  be  treated  statistically. 
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Moreover,  if  SGS  turbulence  is  the  key  for  a  proper  LES 
computation,  the  key  for  an  acceptable  SGS  model  is  its  ability 
to  predict  a  return  of  energy  from  the  subgrid  scales  to  the 
large  eddies.  This  return  of  energy  would  appear  to  necessitate 
new  "de-averaging"  models  when  the  predictive  calculations 
Involve  explicit  computation  for  LES  and  modeling  for  SGS 
turbulence.  This  aspect  could  become  crucial  when  a  vortex  which 
is  subgrid  somewhere  (e.g.,  sufficiently  far  downstream  from  the 
jet  origin)  needs  to  be  explicitly  revolved  elsewhere  (e.g., 
where  the  grid  is  finer),  since  this  implies  the  re-emergence  of 
a  structure  from  the  SGS  model  into  the  resolved  region  of  LES. 

2 . 2  ACCOMPLISHMENTS 


The  different  aspects  of  the  large-  and  small-scale  motions 
addressed  by  the  present  research  are  briefly  discussed  next. 

2.2.1  SGS  Modeling  Procedure 

Current  effort  is  focused  on  one-point  closure  models.  It 
is  worth  noting  that  most  of  the  closure  models  attempted  to  date 
are  either  one-point  or  two-point  models  (depending  on  the  number 
of  spatial  points  appearing  in  the  desired  statistical  results). 
Although  the  two-point  closure  models  are  much  more  complicated 
and  have  been  limited  so  far  to  homogeneous  (and  usually 
isotropic)  flows,  their  application  to  the  circular- Jet 
turbulence  is  clearly  warranted.  Our  initial  effort,  however, 
has  considered  one-point  closure  models.  The  main  motivation  for 
considering  the  candidate  model  is  to  facilitate  the 
incorporation  of  a  variable-viscosity  capability  in  ALGAE  and  to 
provide  for  a  basis  of  benchmark  comparison  with  a  more 
sophisticated  one-point  model  and  possible  two-point  models  that 
will  be  considered  subsequently  in  our  research. 
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As  noted  earlier,  the  averaging  over  the  eubgrld  scales 
introduces  a  coefficient  of  viscosity  for  the  eddy  motion.  Thus, 
the  effective  viscosity  p  Involved  in  the  LES  computations  could 
be  considered  as  the  sum  of  the  molecular  viscosity  of  the  fluid 
n  and  a  turbulent  eddy  viscosity  The  candidate  model  for 

determining  pt  is  the  algebraic  mixing-length  model  due  to 
Launder  et  al.  [1972].  The  formulation  discussed  here  applies  to 
axlsyametrlc  free-shear  flows  with  monotonic  velocity  profiles. 
The  velocity  is  the  maximum  at  the  centerline  and  decreases 
monotonlcally  in  the  radially  outward  direction.  Thus,  no 
recirculating  flow  with  an  inflection  point  in  the  velocity 
profile  can  be  considered.  Clearly,  this  model  is  appropriate 
for  addressing  the  single  round  free  jet.  This  mixing-length 
model  is  one  of  the  three  turbulence  models  implemented  by  Cline 
[1981]  in  the  VNAP2  procedure  for  solving  the  two-dimensional, 
time-dependent,  compressible,  turbulent  flow. 

The  eddy  viscosity,  according  to  this  model,  is  given  by 

Pt  -  Pi2  [  (aw/ar)2  +  (au/az)2] ,1/2 

where  u  and  w  are  the  radial  and  axial  velocity  components,  r  and 
z  are  the  radial  and  axial  coordinates,  and  l  is  the  mixing 
length.  For  monotonic  velocity  profiles,  mixing  length  is 
expressed  as 

I  -  Clrj-rj | , 

where  C  is  a  constant  ( ■  0.11  for  axlsymmetrlc  flows),  and  r2  and 
r1  are  radii  that  correspond  to  certain  arbitrary  values  of 
normalized  axial  velocity.  Thus, 

and 


for  *  °'1- 


for  <«-"alnl/(«.ax-w.ln)  -  0.9 


[10) 


Her*  w  . _  and  w _  denote  the  minimum  and  maximum  axial  velocity 

min  max 

component*  of  the  monotonically  decreasing  velocity  profile.  For 

the  round  Jet  of  our  Interest,  w _ „  is  the  centerline  velocity  at 

max 

a  particular  axial  (z)  location.  For  the  jet  issuing  into 

quiescent  ambient,  wB±n  is  zero.  For  computational  expediency, 

however,  a  nonzero  value  of  0.01  w _ _  (l.e.,  1.0k  of  the 

max 

centerline  velocity)  is  appropriate. 

With  the  typical  boundary  conditions  on  the  centerline 
(viz.,  u»o  and  aw/dr-o,  Invoking  symmetry),  the  equation  for  pt 
must  be  modified  on  the  centerline.  This  is  done  by  replacing 
the  previous  equation  with 

Pt  ■  p*3 |a2«/ar2| 

on  the  centerline. 

2.2.2  ALGAE  Modifications 

The  initial  effort  involving  LES  computations  with  ALGAE  has 
focused  on  certain  modifications  aimed  at  improvement  in  physical 
modeling  and  code  enhancements .  The  former  aspect  Includes  the 
Incorporation  of  a  variable  eddy  viscosity  and  the  refinement  of 
sp  itial  resolution.  The  latter  aspect  Includes  vectorlzatlon  and 
graphic  display  enhancements. 

The  variable  eddy  viscosity  in  ALGAE  is  determined  from  the 
algebraic  mixing-length  turbulence  model  discussed  in  the 
previous  paragraph.  Test  computations  have  been  in  progress  at 
ICMA  to  discern  the  Influence  of  the  variable  eddy  viscosity  on 
the  flowfleld  predictions. 
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The  earlier  version  of  ALGAE  had  a  maximum  number  of  4800 


finite-difference  grid  points  for  the  two-dimensional 
computational  domain.  This  has  been  increased  to  6400  and  is 
expected  to  result  in  more  refined  simulations  of  the  large-scale 
motions. 

The  ongoing  vectorization  involves  the  conversion  of  ALGAE 
(originally  written  for  DEC  1099)  for  operation  on  the  FPS-164  at 
ICMA  and  the  CRAY-XMP  at  the  NSF  Pittsburgh  Supercomputer  Center. 
Further  optimizations  for  the  CRAY  are  in  progress. 

Predictive  calculations  of  the  circular-jet  turbulence  with 
LES  computations  and  SGS  turbulence  modeling  can  be  greatly 
facilitated  by  the  availability  of  enhanced  graphic  display 
capability.  An  appealing  option  is  the  computer-generated  movie 
of  the  evolving  flowfleld.  The  software  necessary  to  generate 
such  16mm  movies  using  the  ICMA  Datagraphlc  Communications  Unit 
has  been  developed  and  is  currently  being  tested. 

2-2.3  Far-Fleld  Boundary  Conditions 

The  specification  of  boundary  conditions  for  the  numerical 
computation  of  the  free-jet  development  presents  many 
difficulties.  Ideally,  our  Interest  is  the  solution  to  the 
problem  of  the  jet  discharging  into  an  unbounded  domain. 

However,  the  present  ALGAE-based  approach  requires  that  the 
problem  be  transformed  into  a  pseudo  boundary-value  problem 
within  the  domain  of  the  computational  grid.  For  Instance,  since 
ALGAE  is  limited  to  treating  confined  flowfields,  it  is  necessary 
to  Impose  far-fleld  conditions  on  certain  artificially  Introduced 
boundaries . 

Initial  numerical  simulations  have  centered  on  a  model  jet 
problem  whose  geometry  is  shown  in  Figure  1.  The  dashed  lines  R 
and  T  represent  the  pseudo  boundaries  in  the  quarter  plane.  Of 
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Figure  1 
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Model  Jet  Problem  Geometry. 
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the  variety  of  conditions  tested  so  far,  specification  of  the 
pressure  on  these  boundaries  appears  to  yield  the  most  realistic 
behavior  of  the  numerical  solution.  Appropriate  pressure 
distributions  may  be  obtained  by  assuming  that  the  far-field 
conditions  are  approximately  those  of  related  model  problems  for 
which  analytical  solutions  are  known.  For  example,  a  linear 
distribution  results  from  a  Poiseuille  flow-type  assumption. 
Another  possibility  currently  under  consideration  (Porsching 
[1986a])  is  to  use  the  pressure  field  corresponding  to 
irrotational  flow  of  a  jet  with  free  stream  lines.  This  is  a 
classic  problem  of  a  type  first  considered  by  Helmholtz  [1868] 
for  the  flow  out  of  a  slot  (plane  jet).  Its  solution  has  been 
reproduced  in  many  hydrodynamics  texts  (see,  for  example,  Milne- 
Thompson  [1955]).  The  solution  technique  uses  a  sequence  of 
Schwarz-Christof fel  transformations  to  obtain  an  implicit 
representation  of  the  complex  potential  of  the  motion. 
Bernoulli's  theorem  may  then  be  employed  to  relate  the  far  field 
pressure  on  streamlines  within  the  jet  to  the  known  pressure  at 
infinity.  Plots  of  the  streamline  and  velocity  fields  resulting 
from  an  application  of  this  technique  to  the  plane-jet  case  may 
be  seen  in  Porsching  [1986a].  Analogous  treatment  of  the 
circular  jet  remains  to  be  carried  out. 

Another  approach  to  obtain  the  far-field  pressure 
distribution  is  through  an  asymptotic  analysis  of  the  fully 
developed  region  of  the  axlsymmetrlc  turbulent  jet.  A  useful 
start  is  provided  by  the  results  of  ongoing  research  of  Bush  and 
Krlshnamurthy  [1986]  which  have  been  shown  to  be  consistent  with 
the  experimental  results  of  Reichardt  [1942]  and  theoretical 
results  of  Gortler  [1942]. 

2.2.4  Weakly  Dissipative  Difference  Methods. 

Sustained  numerical  simulation  of  the  vortex  structures  in 
the  flowfield  of  a  circular  jet  requires  the  use  of  difference 
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methods  with  small  numerical  dissipation.  Unfortunately,  in  the 
absence  (or  near  absence)  of  naturally  occurring  dlaeipatlve 
mechanisms,  it  is  precisely  the  numerical  vlscoeity  that 
stabilizes  the  method.  Thus,  it  is  important  to  examine  the 
prospect  of  constructing  methods  that  are  robust  (with  respect  to 
their  stability  characteristics),  but  not  uniformly  overly 
dissipative. 

A  conceptually  simple  way  to  attenuate  the  dissipative 
effects  of  a  difference  method  is  to  create  a  hybrid  method  by 
blending  the  given  method  with  a  less  robust  but  more  accurate 
one.  The  idea  is  to  design  the  weights  used  in  the  blending 
process  so  that  in  regions  where  little  numerical  dissipation  is 
needed,  the  accurate  method  is  dominant,  whereas  in  regions 
requiring  significant  numerical  dissipation  to  preserve  stability 
(or  certain  qualitative  featuree  of  the  solution  such  as 
monotonlclty) ,  the  original  scheme  prevails.  Such  self-adjusting 
methods  in  computational  fluid  dynamics  were  apparently  first 
considered  by  Harten  and  Zwas  [1972],  and  a  comprehensive  account 
of  the  steps  involved  in  the  deeign  of  these  methods  is  contained 
in  the  thesis  of  Wilders  [1983]. 

Hybridization  is  also  ths  notion  behind  the  Flux  Corrected 
Transport  (FCT)  schemes  of  Boris  and  Book  [1973,  1976]  and  Book, 
Boris,  and  Hain  [1979].  The  FCT  schemes,  originally  developed 
for  one  space  dimension,  were  given  a  nontrivial  multidimensional 
generalization  by  Zalasek  [1979].  He  also  showed  that  they  could 
be  interpreted  in  terms  of  convex  combination  of  flux  terms 
related  to  low-order  (strongly  dissipative)  and  high-order 
(marginally  dissipative)  difference  methods.  The  FCT  weight 
selection  process  uses  a  "monotonlclty  constraint"  on  the 
numerical  solution  and  has  a  particularly  simple  formulation  in 
one  dimension.  Specifically,  the  weights  are  determined  so  as  to 
maximize  the  effect  of  the  high-order  method's  flux  terms, 
subject  to  the  condition  that  over  any  tlmeetep  no  extrema  are 
introduced  that  would  not  also  be  present  in  the  low-order 
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solution  at  the  new  tine.  This  implies  that  the  total  variation 
of  the  hybrid  grid  function  does  not  exceed  that  of  the  low-order 
grid  function.  The  same  Idea  has  been  used  to  define  total- 
varlation-diminlshlng  ( TVD )  schemes  (Harten  [1983],  [1984]). 

The  monotonlclty  constraint  is  consistent  with  the  behavior 
of  a  solution  of  the  one-dimensional  scalar  convection  equation 
[the  total  variation  of  such  a  solution  does  not  increase  in 
time),  and  in  this  case  the  FCT  algorithms  perform  impressively. 
However,  the  situation  is  different  for  systems  of  nonlinear 
conservation  laws.  According  to  Woodward  and  Colella  [1984]: 

Then  no  such  monotonlclty  constraint  is  implied  by 
the  differential  equations,  and  the  use  of  such  a 
constraint  can  lead  to  difficulties.  In 
particular,  a  smooth  region  with  strong  gradients 
can  be  turned  into  a  sequence  of  discontinuous 
jumps,  with  the  appearance  of  a  staircase. 

The  unsuitability  of  a  TVD  condition  in  the  design  of  difference 
methods  for  multidimensional  qumsllinear  systems  is  also 
indicated  by  a  recent  result  of  Rauch  [1986].  There  it  is  shown 
that  unless  the  commutators  of  all  of  the  Jacobian  matrices 
appearing  in  the  system  vanish,  no  multiple  of  the  W1'1  seminorm 
of  the  initial  condition  can  bound  this  seminorm  at  a  later  time. 
Thus,  it  is  unlikely  that  a  numerical  solution  with  a  TVD 
property  will  converge  to  a  solution  of  the  original  system. 

In  view  of  these  difficulties  with  the  TVD  condition,  it  is 
appropriate  to  consider  other  weight  selection  criteria.  One 
such  alternative  is  based  on  the  ability  of  the  hybrid  method  to 
conserve  (or  nearly  conserve)  the  discrete  energy  of  the 
numerical  solution  that  it  produces.  For  the  simple  (constant 
coefficient)  convection  equation 
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In  which  the  stable  upwind-difference  method  Is  blended  with  the 
unstable  centered-dlf ference  method,  this  leads  to  a  hybrid 
method  with  the  single  weight  8  given  by  the  formula 

9  -  [2*  -  1  +  (1  -  4*2)1/2  ]/2A , 
where  A  Is  the  Courant  number. 


When  applied  to  a  standard  square  wave  test  problem,  the 
resulting  scheme  produces  a  numerical  solution  that  virtually 
coincides  with  that  of  the  Lax-Wendrof f  and  leapfrog  scheme. 

This  hybrid  method  is  not  TVD ,  but  does  reduce  the  numerical 
dissipation  of  the  upwind  method  by  about  21%.  Furthermore,  for 
this  problem  it  appears  possible  to  modify  the  analysis  to  permit 
weights  that  dapend  on  the  local  character  of  the  numerical 
solution  (Porsching  [1986b]).  In  this  way  TVD-like  properties 
might  also  be  Incorporated  Into  the  method. 
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